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ABSTRACT 

We use radiative transfer to study the growth of ionized regions around the brightest, 
z=8 quasars in a large cosmological hydrodynamic simulation that includes black hole 
growth and feedback (the MassiveBlack simulation). We find that in the presence of 
the quasar s the comoving HII bubble radii reach 10 Mpc/h after 20 My while with the 
stellar component alone the HII bubbles are smaller by at least an order of magnitude. 
Our calculations show that several features are not captured within an analytical 
growth model of stromgren spheres. The X-ray photons from hard quasar spectra 
drive a smooth transition from fully neutral to partially neutral in the ionization front. 
However the transition from partially neutral to fully ionized is significantly more 
complex. We measure the distance to the edge of bubbles as a function of angle and 
use the standard deviation of these distances as a diagnostic of the isotropy of ionized 
regions. We find that the overlapping of nearby ionized regions from clustered halos not 
only increases the anisotropy, but also is the main mechanism which allows the outer 
radius to grow. We therefore predict that quasar ionized bubbles at this early stage 
in the reionization process should be both significantly larger and more irregularly 
shaped than bubbles around star forming galaxies. Before the star formation rate 
increases and the Universe fully reionizes, quasar bubbles will form the most striking 
and recognizable features in 21cm maps. 
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1 INTRODUCTION 

The current consensus is that the contribution to the global 
ionizing budget from quasars during the Epoch of Reioniza- 
tion (EGR) is small compared to that from early stars and 
galaxies(see, eg, Loeb||2009 Giroux k, Shapiro||1996 



Trac 



Gnediii||2011| . The characteristic proper radius of ioniz- 



ing bubbles at the end of EoR is constrained to be on the 
order of 10 Mpc/h (Wyithe & Loeb 2004), and the quasar 
contribution is limited to be less than 14% (Srbinovsky 
Wyithe 2007| of the total. Even though the contribution 
to global reionization by quasars is constrained in this way, 
bright quasars may still leave a signature on the growth 
of individual ionized regions. Several authors have investi- 
gated this signature in mock 21cm emission spectra taken 
from simulations of isolated quasars ( Majumdar et al.|2011 



Datta et al.|2008| |Geil &: Wyithe|2008t. In observat ions, an 
object recently reported by Mortlock et al. (20TT|, ULAS 
J1120+0641, has a luminosity of 6.3 x W'^Lq ai z ^ 7 and 
a proper near-zone radius of less than 2 Mpc/h. The near- 
zone radius is consistent with the possibility of bright quasar 
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driven growth in a near neutral IGM background ( [Bolton 
et al.||2Qllt . 



In this paper we study the growth of the ionization front 
of bright quasars in an almost neutral cosmological context. 
The quasars and their surrounding medium are selected from 
a large hydrodynamic simulation (the MassiveBlack simu- 
lation, introduced in Di Matteo et al. (2012)), and then 
post-processed with a radiative transfer code. This allows 
us to simulate 8 rare quasars using reasonable computing 
resources. Gur focus is on the evolution and properties of 
the largest individual ionized bubbles, the sources that pro- 
duce them, and the relationship between the two. Because 
the simulation forms quasars and star forming galaxies ab 
initio, we are able to make use of the luminosities and posi- 
tions of the radiation sources that the simulation produces, 
rather than setting them in by hand. In this paper, however 
we do not deal with the full reionization of the volume of 
the simulation, which would require following the evolution 
of the entire density and ionization field from high redshifts 
down to at least z = 6. Instead we restrict ourselves to 
the growth of ionized regions around an early period in this 
process (at ^ = 8), where the photon path lengths are still 
smaller than the computational subvolumes we analyze. We 
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leave the study of the full reionization of the volume to fu- 
ture work. 



2 HYDRODYNAMIC SIMULATION 

The SPH output we use is from the MassiveBlack simu- 



lation (see Pi Matteo et al. (2012); DeGraf et al. (2012); 



Khandai et al. ( 2012 ) for further details), which was run with 
a ACDM cosmology with parameters (Qa, Qm, ^b, Z?^, erg) = 
(0.74, 0.26, 0.044, 0.72, 0.8). A total number of 2 x 3200^ gas 
and dark matter particles were followed in a box of 0.75 Gpc^ 
from red-shift z = 159 to red-shift z = 4.75. This simulation 
is by far the largest cosmological hydrodynamics simulation 
run with the P-GADGET 3 program. 

This run not only contains gravity and hydrodynamics 
but also all the extra physics (subgrid modeling) for star 
formation (Springel & Hernquist||2003 ), black holes and as- 



sociated feedback processes. 

The basics aspects of our black hole accretion and feed- 
back model (Di Mat teo^et al.||2008 ) consist of representing 
black holes by collisionless particles that grow in mass (from 
an initial seed black hole) by accreting gas in their environ- 
ments. The accretion rate is given by 



M 



and 



iin(MB 



MBondi 



(ci+ i;2)3/2' 



(1) 



(2) 



where p and Cs are the density and sound speed of the ISM 
gas respectively, and v is the velocity of the black hole rel- 
ative to the gas, and MEdd = ^Edd/(^c^) (where LEdd is 
the BH Eddington Luminosity, LEdd = 1-26 x 10^^ergs~\ 
At the high redshift and high halo masses we are carrying 
out the analysis here, black holes are accreting at their Ed- 
dington rates ( DeGraf et al.|2012 ) making any detail of the 
subgrid models for black holes virtually irrelevant. 

For the black hole feedback, a fraction of the radiative 
energy released by the accreted material is assumed to cou- 
ple thermally to nearby gas and influence its motion and 
thermodynamic state (typically referred as BH feedback). 
The radiated luminosity, Lr, from the black hole is related 
to accretion rate, Mbh, as 



i^boi - r] (Mbh x 



(3) 



where we take the standard mean value 77 = 0.1. Some cou- 
pling between the liberated luminosity and the surrounding 
gas is expected: in the simulation 5% of the luminosity is 
(isotropically) deposited as thermal energy in the local black 
hole kernel, providing some form of feedback energy. With 
our simulations the activity of quasars is directly derived 
from the accretion history. The luminosity of the stars and 
galaxies can be derived from the star formation rate history 
which is provided by the multiphase star formation model in 
the simulation ( Springel Hernquist|2003| . The multiphase 
star formation model also provides a mechanism to remove 
the self-shielded interstellar medium (ISM) from the matter 
density field. 
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75.5 


34.7 


2.06 


3.93 


1.39 


0.45 


3 





84.9 


7.3 


8.35 


3.71 


0.49 


0.22 


4 


3 


60.1 


119.0 


8.18 


2.72 


7.76 


0.54 


5 





67.1 


16.0 


8.35 


3.71 


0.56 


0.21 


6 


4 


70.8 


5.5 


0.50 


2.72 


0.06 


0.16 


7 


5 


65.6 


13.2 


0.95 


1.98 


0.69 


0.30 


8 


6 


66.6 


6.4 


1.29 


3.82 


0.20 


0.21 


9 


7 


64.9 


7.2 


0.38 


2.22 


0.15 


0.15 



Table 1. UV Flux of Subvolumes. The first column #H (0 - 9) 
identifies the ten most massive halos. The second column #V 
(0 - 7) identifies the eight unique subvolumes that contain the 
halos. The halo mass 

(Mhalo) is in units of 10^0 Moh-^. The 
blackhole mass (Mbh) is in units of 10^ MQh"-*^. The luminosities 
are in units of 10^^ see" Halos 0, 3, and 5 are clustered and 
belong to the same subvolume 0. 
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QUASARS IN MASSIVEBLACK 



The Epoch of Reionization (EoR) is currently constrained 
to be between z = 10, by the CMB measurement (Komatsu 
et al.l 2009), and z = 6, by the Sloan Digital Sky Survey 



measurements of the intergalactic medium surrounding high 
redshift quasars ( |Fan et al.|20Q6l|Willott et al.|2007| . During 
the MassiveBlack run, a uniform UV background radiation 
field is introduced at the end of the EoR in the optically 



thin approximation (see e.g. Bolton & Haehnelt (2007)). At 
z = 8 the uniform radiation has not yet been switched on, 
so that IGM surrounding quasars is mostly fully neutral. In 
a postprocessing step (see below), we then model the radia- 
tion field with radiative transfer, using both a full spectrum 
quasar component for black hole sources and a UV spectrum 
stellar component for the stellar sources. 

In MassiveBlack, the quasars are rapidly accreting su- 
permassive black holes, and their fuelling is driven from the 
large scales and occurs through high density cold flows along 
the cosmic filaments D i Matteo et ah] (2012). In general 
larger halos host brighter quasars unless (i) the quasar is 
turned off by feedback, or (ii) the quasar has not yet grown 
its black hole mass significantly. We select a quasar system 
based on halo mass, taking the 10 most massive halos in the 
simulation at redshift z = S (see Table [3| . For each of the 
halos, we cut a cubical subvolume of side length 50Mpc/h 
from the original simulation. Due to their clustering, three 
of the halos (halo number 0, 3 and 5 in Table [3| form at 
spatial locations within subvolume 0. Only 8 unique sub- 
volumes are identified, which we refer as subvolume to 7 
throughout the rest of the paper. 

Next, we compute the luminosity due to the quasars 
and star formation occuring in each of the subvolumes. The 
assignments of luminosities are described in sections [3T] and 
|3.2| We show the total UV flux of all sources in the region 
and that of the central sources within the IMpc/h radius 
of the center in Table [l] and Figure [l] There is an expected 
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Figure 1. Sample Luminosity. The left panel compares the total UV photon luminosities of the stars and the quasars within the 
subvolumes. The right panel compares the UV photon luminosities from the sources within IMpc/h of the center. The central sources 
directly contribute to the growth of the central ionized region (CIR). The dotted lines represent the linear regression of the data. In the 
bottom panel the central luminosity as a fraction of the total luminosity is shown. The subvolumes fall into three catogories: (S) Stellar 
Driven, (QS) Subdominant Quasar, (Q) Quasar Driven. 



correlation between the central quasar and central stellar 
luminosity. Based on the ratio of the central luminosity to 
total ionizing luminosity within the subvolume, as well as the 
composition (quasar or stellar) of the central luminosity, we 
can divide the subvolumes into three groups: 

• Stellar (S): The quasar luminosity in the central sources 
is smaller than the sum of stellar luminosities, and the entire 
central source luminosity is small comparing to the ionizing 
luminosity of the entire box. We note there can be two pos- 
sible distributions of sources within the box: one being the 
presence of another major source that is not near the cen- 
ter of the subvolume; the other being ionizing sources that 
are rather more uniformly distributed throughout the sub- 
volume. In our subvolumes, the latter is more often to the 
case. 

• Subdominant Quasar(QS) : The central quasar lumi- 
nosity dominates the central luminosity, but is not a sig- 
nificant fraction of the total luminosity of sources in the 
subvolume. 

• Quasar (Q) : The central quasar luminosity dominates 
the entire subvolume. 



3.1 Quasar Radiation 

We calculate the bolometric luminosity Lboi of each quasar 
from the accretion rate M of the supermassive black hole 
with equation (pi. We choose to use as M a single value 



Band Range 



UV 


13.6 eV to 250 eV 


Soft X-Ray 


250 eV to 2KeV 


Hard X-Ray 


2 eV to 10 KeV 



Table 2. Band definitions. 



throughout the calculation, which is averaged accretion rate 
measured from the simulation over a timescale tq — 2 x 
lO^yrs, where tq is also the length of time over which we 
follow the evolution of the quasar ionization front. We esti- 
mate the luminosities in different bands from the bolometric 
luminosity according to the fitting formula of [Hopkins et al.] 
(|2007|). The band definitions are listed in Table |2| The lu- 
minosity in UV band is estimated from B band luminosity, 
assuming power laws with optical index of auv = 1-76 and 
as = 0.44 respectively, 

/ 250 \ 1~"UV 

i„„ = ^ C^)-^ (EL)-"^ %) Zl, (4) 

where ux = c/120nm is the pivot from optical to UV, ub = 
c/445nm and lyj — c/91.1 nm. We calculate the photon flux 
of the quasar sources according to these power law spectra 
applied to the various bands. 
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3.2 Stellar Radiation 



5 RADIATIVE TRANSFER SIMULATION 



The stellar luminosity of a halo is given by 

iVsFR = fescN^/bXH^/rUp, (5) 

where $ is the star formation rate of the halo, Xh = 0.76, 
and rrip is the proton mass. We take the escape fraction 
to be /esc = 0.1, and the photon to baryon ratio to be 
N^/b = 4000, a value found to produce a reionization his- 



For this study we have rewritten the SPHRAY Monte Carlo 



tory consistent with observations by Furlanetto et al. ( 2004 ) ; 
[Sokasian et al. (2003 ) . For simplicity, the stellar sources are 
assigned to the center of the halo; in future work we plan to 
investigate the effect of the positioning of the stelllar sources 
by either following the star forming gas particles in the sim- 
ulation directly or the gravitationally bound subhalos. To 
ease the modelling, we only include and assume the same 
UV band spectrum for the stellar radiation to that of the 
Quasar spectrum. 



4 ANALYTIC GROWTH OF STROMGREN 
SPHERES 

The growth of Stromgren spheres in a clumpy cosmologi- 



cal environment can be analytically modeled using (Cen & 
iHaiman 2000), 



3N, 



ph 



(6) 



dt ' 47r(nH) 

where the first term represents the Hubble expansion, and 
can be neglected in the current context. The solution is 
straightforward (eg, Majumdar et al.]|2011 ) 

-1/3 / ^ \ V3 



Rt 



Rs 



(l - exp (- 



ih) 



(Ch (nn) as) 
3iVph 



(7) 



1/3 



We use this analytic solution to compare to our numeri- 
cal results, but first we need to measure the constants in the 
above expression from the simulation sub volumes. We cal- 
culate the clumping factor Ch using a method appropriate 



for an SPH simulation, described by Pawlik et al. (2009), 



Ch 



(8) 



where hi is the SPH smoothing length. The clumping factor 
we measure from the gas in the entire subvolumes has a 
value Ch ^ 4, giving a recombination time about ts ^ 
lOOmyear ^ tq. However we note the clumping factor can 
go up to ^ 20 within IMpc/h radius from the center of 
the subvolume, reducing ts to ~ 20myear. For the values of 
pi, we use the density of the intergalatic medium(IGM) gas 
component which is relevant to the radiative transfer: the 
interstellar medium (ISM) is excluded in consistency with 
the simulation; see below. 

For the source luminosity, we make the simplest ap- 
proximation, summing the luminosities of all sources with 
in the subvolume to obtain one effective iVph. We also as- 
sume that the volume consists of pure hydrogen (Xh = 1, 
i.e. no helium) in the calculation of the analytic model. The 



radiative transfer code|Altay et al. 2008 Croft & Altay 



2008) in C with OpenMP(TM), so that it runs in parallel 



on shared memory systems. We have also eliminated the on- 
the-spot approximation of recombination, so that instead re- 
combination rays are emitted when the recombination pho- 
ton deposit at a spot exceeds a threshold, as done in the 



code CRASH(Masem et al. 2003). The parallel version al- 



lows us to trace many rays at the same time on several CPUs 
within a single time step, achieving a near- linear speed up 
with a small number (10~30) of CPUs. Monochromatic rays 
are sampled from frequency space according to the source 
spectrum. Hydrogen and Helium ionization and recombina- 
tion are both traced. The tabulated atomic reaction rates of 



Hui & Gnedin ( 1997) taken from the serial version SPHRAY 



are used. 

As a new addition to the code we model the secondary 
ionization of high energy photons using a fit to the Monte 



Carlo simulation results of Furlanetto & Stoever (2010) 



Shull & van Steenberg (1985). These fast non-thermal 



recombination rate as is taken from |Hui Gnedin (1997) 



equilibrium electrons produced from high energy ionization 
photons (e.g.. X-ray photons) may collide and secondarily 
ionize more neutral atoms, increasing the ionizing efficiency 
of the quasar sources by allowing one high energy photon 
to ionize tens of neutral atoms. On the other hand, our cur- 
rent post processing approach is incapable of modeling the 
heating from the residual energy of the ionizing photons, 
because in this approach the thermal evolution is decoupled 
from radiative transfer. MassiveBlack contains a calibrated 
thermal term in the black hole feedback mechanism, which 
addresses this issue, at the expense of injecting the heating 
energy only locally to the vicinity of the black hole. 

The inter-stellar medium (ISM) component of gas par- 
ticles condenses into cold clouds which are self-shielded from 
cosmic radiative transfer (except that they may host stellar 
radative sources) . As a result we excise them from the mat- 
ter density field when performing the ray tracing calculation 
(by doing this we also assume that their small cross-section 
would not have affected ray tracing for the rest of the gas). 
MassiveBlack has a multiphase (cold, hot gas phases) star 
formation model, from which we identify the cold cloud com- 
ponent with the ISM (S pringel &: Hernquist|2003t . The cold 
fraction is an increasing function of the mean density of the 
gas particle, effectively removing the densest particles from 
the density field in a manner similar to the threshold method 
used to calculate the clumping factor by |Pawlik et al.| ( |2009| ), 
and, the removal of the cold high density gas in the X-ray 
emission calculation by [Croft et al.| ( |2001| ) . 

In MassiveBlack the total(IGM + ISM) mean baryon 
density around quasar sources can be as high as 60 cm""^. 
Were the ISM not removed from the ray tracing, the high 
mean density would shield off the radiation and prevent the 
growth of any cosmic scale ionized regions. Figure [2] shows 
the histogram of the overdensity of gas particles before and 
after the removal of cold clouds in one of the subvolumes. 



5.1 Parameterizing the Ionization Front 

We use the center of the chosen halo to define the center 
of the subvolume, and also the central ionized region (CIR). 
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Figure 2. IGM density and total gas density histogram The 
dashed Hne is the total gas density, while the dotted line is the 
IGM density. The histogram is performed on the overdensity of 
the particles in the simulation. Using the IGM density results a 
reduction in the number of dense particles and allows a sensible 
treatment of the radiation transfer. 



When quantifying the size of the CIR we measure the spher- 
ically averaged fraction in a species, x using (we show the 
example of HII regions, similar definitions hold for other 
species): 



XHi(r) 



J6{r' -r)xH\pdV ^ E.^en 
J6{r' -r)pdV T.sheii'^ 



(9) 



Here the sums are carried out by binning SPH particles (in- 
dex i) in shells. Particles are assigned as a whole to shells by 
their center position rather than being intergrated over the 
SPH kernel with in the shell. We define the averaged CIR 
radius at a given threshold x* to be the first crossing of x(r) 
with X*, or 



Rs{x) 



if{r|xHi(r) = X*}. 



(10) 



The first crossing corresponds to the edge of the CIR, and 
later crossings indicates the edges of the nearby ionized re- 
gions. The motivation for this definition comes from its simi- 
larity to that used to quantify ionized regions seen in quasar 
absorption li ne spectra, the Lv o^ absorption Near-Zone ra- 
dius Rnz by Fan et al. (2006). In practice, we take as Rs 



the average of the two nearest bins tl and tr which satisfy 
[xH\{rL) - x\[xH\(rR) - x] < 0. 

The spherically averaged radius Rs is a simple quan- 
tity to use in comparisons but the angular variations in the 
properties of the CIR are completely lost. To quantify the 
angular dependency, we bin the SPH particles into angular 
cones and calculate the averaged fraction within the cones. 



XHi(r, 6',(/)) 



y 



, XH\,imi 



E 



cone. shell 



The angular dependent CIR radius is 



Rs{x,e, 



inf{r|xHi(r, 6*, 



(11) 



(12) 



Note that in general < Rs{x,0,(j)) >^ Rs, because averag- 
ing and binning do not commute. For measuring i?s(x, 0, 0), 



we use 12 x 16^ cones. We note that increasing to 12 x 32^ 
cones does not significantly alter the results. 

We measure the anisotropy of the ionized bubble from 
the variance of Rs 



As{x) = STk[Rs{x,Ok 



(13) 



where STk stands for the standard derivation. We run the 
simulation with sufficient number of rays so that the typical 
shot noise contribution towards As is negligible, as shown in 
the next section. 

Three different levels of species fraction x are used to 
define three different levels of CIR fronts in this study: 

i Inner front, x = 0.1 for the neutral fraction, or x = 0.9 
for the ionized fraction; 

ii Middle front, x = 0.5 for the neutral fraction; 

iii Outer front, x = 0.9 for the neutral fraction, or x = 0.1 
for the ionized fraction. 

The Inner front corresponds to the near-ionized edge of 
the Stromgren sphere, and the Outer front correcponds to 
the near-neutral edge of the Stromgren sphere. We however 
note that these choices are for illustrative purposes only and 
they do not have any direct correspondence with threshold 
values used to detect 21 cm or Lya observation signatures. 



5.2 Shot-noise and Convergence 

One issue which must be addressed in Monte Carlo ray trac- 
ing (RT) schemes is the prescence of shot noise and its ef- 
fect on convergence of results, something particularly im- 
portant when sampling is also in frequency space (eg ,Iliev| 
et al.|2006| on CRASH). Shot noise artificially increases the 
angular anisotropy As measure, and so is of direct concern 
for this study. 

We define a shot noise parameter 7 to be the ratio be- 
tween the number of photons in a ray packet, n° and the 
number of atoms in an SpH particle n^. 



(14) 



A small 7 guarantees that the ionization front cannot ad- 
vance by more than one particle in one time step. When 
7^1, the ionization front advances slowly and the shot 
noise is controlled. One of course still needs to ensure the 
rays have a sufficient angular resolution to resolve the angu- 
lar scale used in the calculation of the anisotropy. 

In MassiveBlack, a typical SpH gas particle has a mass 
of 5 X lO^Mo/h, equivalent to = 9 x 10^^. For ray trac- 
ing with 10^ steps and 128 rays per time step, (1.28 x 10^ 
rays), and with the total luminosity listed in Table [l] an av- 
erage packet contains n° = 3 X 10^^ photons. The shot noise 
parameter is therefore 7 = 0.03 ^ 1 for our typical runs. 

We can also empirically confirm that convergence is 
reached by increasing the number of rays used in the simu- 
lation and showing that the quantity of interest is insensi- 
tive to further increases. We perform such a convergence test 
with 3 runs on subvolume 4 with (i) 1.3 x 10^ rays, 7 = 0.003, 
(ii) 1.3 X 10'^ rays, 7 = 0.03, and (iii) 1.3 x 10^ rays, 7 = 0.3. 
We quantify the similarity of the ionization front using the 
correlation coefficient of RsiO^cj)) between runs, shown in 
Table [3] The higher correlation between the runs with more 



6 Yu Feng et al 



Runs 


HII Inner 


HII Middle 


HII Outer 


0.003 and 0.03 


0.96 


0.98 


0.86 


0.03 and 0.3 


0.74 


0.94 


0.68 



Table 3. Correlation Coefficients. Shown in the table are the cor- 
relation coefficient of the center bubble radius Rs(0,(f)) between 
different runs. 
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Figure 3. CIR Growth in a Uniform Density Field. The symbols 
are: (i) cross[+]: UV only; (ii) cross[x]: UV and X-ray (iii) circle: 
UV and X-ray, with Secondary ionization. The analytic model 
(dashed lines) neglects Hellium and assumes Xh = 1.0. 



streaming power law with an index of 1/3 in equation ([7|. 
The inner and middle fronts are smaller than the analytic 
model due to a lower ionization threshold used to define the 
front, (ii) The effect of X-ray photon and secondary ioniza- 
tions is quite limited on the inner and middle front. The 
effect is more evident on growth of the outer front. This is 
because the secondary ionizations affect the way the harder 
photons interact with the matter the most, and harder pho- 
tons travel further into the neutral region. We note that by 
adding in the effect of secondary ionizations increases the 
outer radius of ionized bubbles by approximately 10%. 



6.2 MassiveBlack: Visualization 

We show the CIR bubble of the most extreme Q type and 
S type subvolumes( sub volume number 3 and 4 in Table [s] 
respectively) in Figure [4] The halo mass of both are simi- 
lar, 60 X lO^^Moh-^ and so are their total star formation 
rates. However subvolume 3 has a large ionized HII region 
associated with the bright active quasar in the center of the 
subvolume, whilest in subvolume 4 the HII region formed 
merely from stellar sources are small. Note that, even though 
of comparable mass halo, the central black hole is much 
smaller in subvolume 3 and hence its activity has a negli- 
gible impact. This example aims to illustrate how, in the 
presence of an active early quasar, the ionized bubble is far 
more extensive than one driven by star formation alone. 

The time evolution of the HII bubble in subvolume 0, 
where three halos coexist, is shown graphically with a slice 
through the center of the simulation volume every 2 x 10^ yr 
in Figure [5] It is interesting to observe that a neighbouring 
bubble on the right merges with the center bubble and this 
increases the size of the outer front in that direction, con- 
tributing to the anisotropy of the front. We investigate this 
issue in the next section in a more quantitative way. 



rays indicates that the simulation has effectively converged. 



6 RESULTS AND DISCUSSION 
6.1 Uniform Density Field 

We first test P-SPHRAY with a source in a uniform density 
field at z = 8, with H number density riH = 2 x 10~^ cm~^, 
and uniform temperature T = 10^ K. The hydrogen mass 
fraction is Ah = 0.76, the cubical box has a side length of 
50Mpc/h, and we evolve the radiation and ionization frac- 
tions for a duration of 2 x 10^ yrs. We carry out 3 separate 
simulations, considering different central source properties 
for each one as follows: (i) 1.2 x lO^Sec"^ UV; (ii) 1.2 x 
10^^ sec-^ UV, 0.3 X 10^^ sec"^ soft X-Ray, 0.5 x 10^^ sec"^ 
hard X-Ray; (iii) same as (ii), with secondary ionizations. 
The luminosity of the source is motivated by the luminosity 
of the center sources in subvolume 0. The growth of the CIR 
fronts in the simulations are shown as a function of time in 
Figure [3] We show as separate panels the behaviours of the 
three parts of the ionization front: inner, middle and outer 



(as defined in section 5.1). 

There are two highlights from the uniform density field 
simulations: (i) The outer front agrees well with the free 



6.3 MassiveBlack: Spherically Averaged Radius of 
Ionized Regions 

We now examine the results of the radiative transfer post- 
processing of the eight subvolumes from MassiveBlack. The 
averaged CIR radius Rs for HII and Helll, as defined in sec- 
tion |5.1| are shown in Figure |6] as functions of the central 
source flux. We also plot fits of the growth to the 1/3 free 
streaming power law ([7|. At the low luminosity end (mainly 
S type and QS type subvolumes) there is a significant devi- 
ation from the fit for the inner and middle fronts. The outer 
front is more extended than the simple uniform field simu- 
lation, although the source spectra are similar, hinting that 
the structure in the IGM is also contributing to the smooth- 
ing of the front. This smoothing is a phenomenon similar to 



that described by Wyithe & Loeb ( 2007 ) . We also attribute 



the effect to the clustering of sources; however unlike the 
smoothing due to secondary ionization, the clustering con- 
tribution is anisotropic, and we discuss it in section [6^ 

In Figure [6] we can see that the HII and Hell CIR radii 
appear to grow together. The correlation between the two 
can be seen directly by plotting one against the other, which 
we do in Figure |7| where we find there is a strong correlation 
between the HII radius and Helll radius. The Helll radius 
is smaller than the HII radius, agreeing with the finding of 
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Figure 4. Q (Quasar, left) type and S (Stellar, right) type ionized bubbles in MassiveBlack. Both panels are of comoving length 15 Mpc/h 
per side. A camera is put at about 60 Mpc/h from the center of the subvolume and a perspective transformation is applied to form a 
projection onto the image plane of the camera. Red color corresponds to fully neutral IGM and blue color corresponds to fully ionized 
IGM. Yellow is in between two states. Crosses are the sources; both stellar and quasar sources of the entire subvolume are shown. It is 
interesting to notice the lack of a major source (comparing to the bright quasar in subvolume 3) in subvolume 4. 
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Figure 5. Evolution of the ionization front. The evolution of the ionization front in subvolume #0 is shown. The ionization front(a;Hl ^ 
[0.1, 0.9]) is marked in black, and the fully ionized near zone(a::Hl < 0-1) is marked with gray. 



Friedrich et al. (2012) who used the ray tracing code C2Ray. 



We note however that the treatment of secondary ionization 
and Helium ionization in P-SPHRAY is similar to that of 
C2Ray, and an agreement is not surprising. 

The time evolution history of the HII CIR radius for all 
subvolumes is shown in Figure [S] We compare the evolution 
of the three parts of the fronts with the prediction of the 
analytic model(in equation [7|, in which we have used the 
clumping factor and mean H density measured within the 
final CIR in the simulation to estimate a fiducial recombi- 
nation time is- 

The growth of the CIR in the Q type subvolumes flat- 
tens off much earlier than one might expect from a free 
streaming law. In order to ascertain the relevant physical 
timescales, we fit the time evolution of the fronts to the full 
analytic model in equation ([7|, and extract the effective re- 
combination time is as a fit parameter, then mark the time 
in the plot. We can see that for the quasar driven (type Q) 
subvolumes, all three fronts (inner, middle and outer) growth 



stop at about lOmyear, on the same order of the life time 
of the quasar tg ~ 20myear and the fiducial recombina- 
tion time is ^ 20 myear. The deviation from free streaming 
indicates that by merely increasing the life span tq of the 
central quasar, we can not substantially increase the size of 
their CIR. 

In S type type subvolumes, there is a similar tailing 
off for the inner fronts. However, the outer fronts continue 
their free streaming growth, behaving differently from the 
analytical model, which has stopped due to reaching the 
recombination time. We attribute this apparent excessive 
growth of the averaged S type outer front to the anisotropic 
growth via merging with other small ionized bubbles that 
are close to the CIR, as described later in the paper. 

We note that our Monte- Carlo RT approach allows for 
superluminal growth of the CIR. This unphysical situation 
happens well before the first snapshot time which is 2 mil- 
lion years after the sources are turned on. The first snapshot, 
in which the stromgren sphere typically has grown to about 
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Figure 8. Evolution of three Quasar driven subvolumes. The top panels show the CIR radius as function of time; the thick black line 
is the fiducial analytic model as described in the text (equation [?]). The vertical short lines show the recombination time in an analytic 
model as if it is fitted as a parameter to the simulation data points, which are marked with squares (Inner), triangles (Middle), and 
circles (Outer). The bottom panels show the ratio between the radius and the fiducial model. 
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Figure 6. HII and Helll CIR Front Radius. The left panel is for 
the HII CIRs, the right panel for the Helll CIRs. The horizontal 
axis is the central source flux. The dotted line in the flrst two 
panels is fltting to the i/s free streaming power law. 



10% of the full size, gives an upper bound to the contribu- 
tion from super luminal growth. We conclude that the su- 
per luminal growth occuring at early time(< 2 million years) 
and small radii(< 10%) does not play an important role in 
the final shape of the ionized bubble. 



6.4 MassiveBlack: Anisotropy of Ionized Regions 

In Figure [5] we show the anisotropy measured using equa- 
tion ( |13| of all 8 subvolumes. The anisotropy of different ion- 
ization fronts as a function of radius are marked with differ- 
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Figure 7. Correlation between Helll CIR and HII CIR. Helll 
CIRs is plotted against HII CIRs, and the linear regression for 
three fronts are shown with the dotted line. 



ent symbols. As a reminder, in equation (13), the anisotropy 



As is the standard deviation of the ionization front radius, so 
that by comparing to the radius plotted on the x-axis it can 
be seen that the fronts are not very far from spherical sym- 
metry (~10% variations). We find that larger ionized regions 
(corresponding to brighter quasar sources) are in general as- 
sociated with more anisotropy. 

In the simulation, the anisotropics of the inner and mid- 
dle radii of the HII regions do not significantly depend on 
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Figure 9. Anisotropy .vs. Radius. The left panel shows the 
anisotropy of three (inner, middle, outer) HII CIR fronts defined 
as in equation ( |13| ). The right panel shows the anisotropy of the 
Helll fronts. The dotted lines are fits against a square-root + 
linear offset model and are only meant to guide the eye. All 8 
subvolumes are displayed. 



the radius (the curve is relatively flat), but the outer fronts 
have more anisotropy. The Helll regions show a similiar fea- 
ture, except for the inner fronts of three type Q subvolumes, 
which have significantly stronger anisotropy. 

These phenomena lead us to the following explanation, 
incorporating two contributing mechanisms: 

(i) The anisotropic distribution of gas that attenuates the 
ionizing photons; when the density is high in a particular 
direction, the extra absorption decreases the radius of the 
ionized region in that direction. 

(ii) The merging of nearby bubbles from clustered halos; 
when the density is sufficient to host bright sources lying 
in a certain direction, their extra photoionization increases 
the bubble radius in that direction. The outer front is more 
sensentive to merging than the inner and middle front. 

For a small CIR(or a CIR in its early growing stage), 
no merging has occurred and only the density induced 
anisotropy is present. As the CIR grows, it overlaps 
with nearby ionized regions, resulting the second type of 
anisotropy (due to overlaping). We note that this does not 
contradict the finding that clustering contributes to the 



smoothness of the extended ionization front, (due to Wyithe 
& Loeb (2007)), because after spherical averaging, the over- 



all effect on the merged bubbles is to produce smoother ion- 
ization fronts. 

In order to better visualize the structures which cause 
these anistropies, we plot the distance to the different part 
of the ionization fronts in a Molleweide projection as seen 
from the point of view of the central source. These maps of 
the angularly dependent HII bubble radius, Rs are shown 
in Figure [To] where we plot the quasar driven subvolume 
and stellar driven subvolume 4. 

Looking from the top panels downwards for subvolume 
0, we can first see clumps of mostly neutral gas close to the 
quasar which restrict the distance to the xhi = 0.1 fronts 
to be very close by (< 1 Mpc/h), compared to the mean 
distance for this neutral fraction of 6 Mpc/h. In the xhi = 
0.9 plot, for subvolume a prominent red region can be seen 



at the top. This corresponds to a bubble which has merged 
with the central bubble. Next to each panel in Figure [lO] we 
show a histogram of the Rs values. The secondary bubble is 
responsible for the small tail of high values region in xhi = 
0.9 histogram, as well as giving an extra contribution to As. 

Moving on to the stellar driven subvolume 4, right hand 
panels in Figure^] we can see that the ionized region radius 
is smaller by approximately a factor of 3. Because this radius 
is small compared to the distance to nearby major halos, 
there is no sign of merging with other large bubbles, and 
the ionized region remains more isotropic. In the bottom 
panels of Figure [lO] we show the Molleweide-projected mass 
density within the inscribed sphere of the subvolume. It is 
interesting to compare this dumpiness with the structure in 
the ionized bubble radius. We can see some correlation with 
some structures, but not as much as might be expected, 
indicating that the interaction of radiation with the clumpy 
medium surrounding the sources is a complex process. 

In Figure we show an orthogonal view to Figure |10[ 
showing structure along rays moving out from the central 
source in some different directions. The outer HII front radii 
in three chosen directions are shown as a function of time. 
This allows us to see in detail how the clustering of halos 
affects the center bubble radius through merging as the sim- 
ulation develops. The three directions are (i) one directed 
towards the second brightest source in the subvolume , (ii) 
one directed towards the merging bubble seen on the right 
hand side of Figure [5] and finally (iii) one oriented in a 
random direction not crossing any nearby secondary bub- 
bles. In order to show the stronger bubble merging effect 
which would happen with more luminous sources, on the 
right hand side of Figure pT] we plot the same 3 sightlines, 
but after increasing the stellar luminosities by a factor of 10. 

Merging happens as the central bubble touches the sur- 
rounding ones and meanwhile the radius along that direc- 
tion significantly increases, compared to the radius without 
a nearby halo. This can be seen most clearly for sightline (iii) 
on the right in Figure [TT] The development of the anisotropy 
due to merging of nearby ionized bubbles is responsible for 
the apparent rapid growth of the averaged outer Stromgren 
radius after the growth of the inner radius is stabilized for 
S type and QS type subvolumes. As the neutral fraction at 
the overlapping edge slowly drops below the threshold, two 
bubbles merge and the averaged radius increases. However, 
at the luminosities seen in the simulation (left hand panels 



of Figure 11 ), such merging mechanism is mostly limited in 
the growth of the outer front, for the slow growth of the 
pre-merging bubbles. It is also worth noting that in the ab- 
sence of the effect of smoothing due to secondary ionization 
and the steep bubble edges will make the merging even more 
difficult to happen. 



7 CONCLUSIONS 

We have presented results from radiative transfer simula- 
tions in the vicinity of high redshift quasars in the Massive- 
Black simulation. We find that the rare brightest quasars 
drive a much more significant growth of ionized regions than 
in the purely stellar driven case. The ionized regions associ- 
ated with active quasars are characterized by (i) a smooth 
ionized fraction transition from the middle to the outer 
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Figure 10. Angular Distribution of the HII Radius. The left column is subvolume 0(Q type) and the right column is subvolume 4(S 
type). From top to the bottom the inner (xm^o.i)? middle (xmrro.s)? outer (xhi=o.9) fronts are shown. The histograms measure the 
variance of Rs in different angles, characterizing the anisotropy. The red region in the xhi = 0-9 plot of subvolume is due to merging 
with a further away HII bubble. The less significant red in the center corresponds to the merging shown in Figure [s] 



front, and (ii) an increased anisotropy in the front when 
it starts to overlap the nearby ionized regions. The nature 
of such growth is significantly more complex than a sim- 
ple analytic growth of a single center bubble with clumping 
correction. 

The largest HII bubble obtained in this simulation has 
a comoving radius of lOMpc/h, which is smaller than the 
general expectation that can fulfill the reionization of the 
universe. ( [TVac Gnediri|[20TT| [Wyithe Loebl[2004l paT 



jumdar et al.|2Qllt The quasar near zones we have presented 
in this paper are the primordial ancestors of the later much 
larger Stromgren spheres which will form near the end of the 
EoR. They are however relatively isolated regions that could 
be interesting objects for study in future 21cm surveys. Af- 
ter the z = S epoch we have modelled in this paper, we 
expect that the global star formation in the MassiveBlack 
simulation increases significantly. This will eventually lead 
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Figure 11. Time evolution of a^Hl- Time evolution of xhi in sub- 
volume 0, along three selected sightlines: (i) the blue line is di- 
rected towads the second brighted halo in the subvolume; (ii) the 
green line is directed towards the merging bubble on the right 
shown in the slice plots; (iii) the red line is directed towards a 
random direction with no nearby bubbles. The left panel is the 
luminosity model used throughout the paper, while on the right 
the stellar luminosity is boosted by a factor of 10 to demostrate 
a stronger merging effect. 



to the global reionization of the universe. We plan to study 
this process and role of quasars in future work. 
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